Journal of Computational Biology

Integrative analysis of the properties of multiple networks pertaining to a biological system is a key problem in systems biology. But computing the properties of thousands of large networks poses a challenge. Most techniques that partially address this challenge, such as parallelism and optimized algorithms, treat the analysis of each network as an independent task. However, the networks arising from the same biological system are often similar to each other, and can we harness this similarity to achieve speedup? Towards this end, we propose a phylogeny-like data structure, nc-tree, to compactly represent an ensemble of related biological networks; and a versatile framework, ProcessNets, that runs incremental algorithms over the nc-tree to speedup property computations on the input networks. Our asymptotic analysis of space and time complexity, coupled with empirical evaluations of diverse ensembles of simulated and real-world GTEx gene coexpression networks, demonstrate that ProcessNets leads to significant space and time advantages when computing different network science measures on ensembles of similar networks. For instance, we achieve a compression factor of 2.1x and 8.1x in storing two representative ensembles of GTEx Muscle Skeletal networks with 1000 nodes/genes, and a time speedup of 2.86x to 3.97x over the fastest baseline for computing degree centrality in ensembles of networks with at least 14,500 nodes derived from subsampled GTEx datasets of ten tissues. These results are promising and encourage application of ProcessNets to analyze biological network ensembles derived from rapidly accumulating consortium/biobank-based datasets.

In network biology, analyzing pathways is critical to study mechanisms or sequences of biological events. The pathway co-membership approach aims to link proteins belonging to the same pathway inside a network that is a clique. Although this method is great for projecting pathway knowledge onto proteins, it remains unsuitable for some downstream applications relying on graph distances like random walks with restart. To address this limitation, we introduce a weighted pathway co-membership method that incorporates biological distances between proteins. We achieve this by developing an abstraction framework for Reactome pathways, which captures both their hierarchical organization and sequential relationships. By assigning weights to these abstractions, we define a scoring system for proteins within the pathway co-membership clique. We illustrate our approach using a Reactome pathway ("Signaling by TGF-beta receptor complex", R-HSA-170834), and explain how we generate the weighted co-membership cliques of the 29 Reactome top-level pathways. Next, we construct the entire weighted pathway co membership network of the Reactome knowledge base. Validation and evaluation of our method is performed using cancer hallmark gene sets and a random walk with restart experiment, comparing our approach to state-of-the-art methods. We show that our approach outperforms state-of-the-art methods for gene prioritization on cancer hallmarks.

Next-generation sequencing data is compositional in nature and ignoring this fact can significantly impact downstream data analyses. This is especially true in co-occurrence (correlation) networks because compositional data is prone to many false negative correlations. The centered log-ration (clr) transformation is a commonly used compositional data analysis (CoDA) approach that converts compositionally constrained data into real number space. The clr transformed values can subsequently be used in machine learning and statistical analyses such as Spearman or Pearson correlations used to generate co-occurrence networks. Here we show that the necessary zero substitution performed prior to clr transformation can create numerous multicollinearity artefacts. We also show that a newer CoDA approach called the robust centered log-ratio (rclr) eliminates these artefacts in correlation-based network analysis.

We have developed a method to integrate transcriptomic coexpression networks across diverse experimental conditions within a single species. Our framework expands previous pannetwork approaches into a hypergraph based pannetwork. It first identifies coexpressed gene clusters within each individual dataset and extracts frequently co-expressed gene sets across multiple datasets using a frequent itemset mining algorithm. This process yields a hypergraph where each hyperedge is assigned a frequency (universality, U). We then extracted a subnetwork comprising high U hyperedges as the core network and formed modules within it. We applied our method to 106 Escherichia coli transcriptome datasets from the GEO database. The modularity of the core network peaked at a universality cutoff of 15, which was subsequently used to define it. Approximately 70% of the resulting core modules correspond to operons, and conversely, approximately 70% of all operons are covered by these core modules. We visualized the core network via an inter-modular network and analyzed core module dataset relationships using a modularity profile matrix. Based on these analyses, we successfully visualized the dynamic reorganization of the coexpression network in response to environmental changes in bacteria.

Single-cell sequencing provides a rich source of data that can be used to model and estimate the evolution of somatic cells within individuals. Many common models for modeling DNA evolution in phylogenetics are based on continuous-time Markov processes that operate on aligned sites. For haploid DNA, such as found in mitochondria, there are many well established models for these processes but since nuclear DNA is diploid, containing copies of both the maternal and paternal genome, it is natural to develop continuous-time Markov processes adapted to diploid DNA. One possible model is that both the maternal and paternal evolve independently according to haploid models so that, in principle, DNA comes with a phase. i.e. a label that tells if the site was derived from the maternal or paternal sequence. However, this information may not be available in practice and only an unphased diploid genotype is available. That is, an unphased model groups both of the possible maternal-paternal labels into a single state. This grouping process, called ``lumping'', is used by other methods as well and results in a coarser process. Care must be taken so that the coarser process can still be modeled as a continuous-time Markov chain; the lumped process of the original phased Markov chain may or may not be Markovian, depending on restrictions imposed on the generator matrix. Here, we explore the mathematical formulations of phased and unphased Markov processes, as well as a collection of other diploid DNA models from the literature.

Learning from imbalanced data remains a major challenge for graph neural networks (GNNs), as minority nodes are not only rare but also structurally marginalized within the graph. We address this issue with CLARA, a hierarchical learning framework that decomposes node classification into two stages: a coarse subgraph-level classifier that selects regions likely to contain minority instances, followed by a fine-grained node-level predictor within these regions. This design improves sensitivity while maintaining scalability, filtering out irrelevant areas and focusing learning on topologically meaningful neighborhoods. Experiments on benchmark graph datasets demonstrate substantial gains over established imbalance-handling methods, with CLARA reaching an F1-score of 88.3\%. The same strategy achieves significant improvements in protein-ligand binding site prediction, underscoring its broad and consistent effectiveness across both biological and general graph learning tasks.

Single cell RNA sequencing (scRNAseq) has revolutionized cellular biology over the past decade. At its most basic, a scRNAseq dataset is a cell by gene matrix with entries giving the transcription count of a particular gene in a particular cell. A common workflow applies log-normalization and z-scoring to the count matrix and then embeds the cells in a lower dimensional space using principal component analysis (PCA). A large number of computational tools use the PCA embedding as a starting point, but statistical models for the distribution of cells in the PCA space are lacking. As a result, many tools depend on ad-hoc assumptions that likely lead to signal loss. In this work, we address this gap using the metacell model, a statistical model of count matrices introduced by previous authors. We show that application of the workflow to count matrices generated under a metacell model yields PCA embeddings distributed as a Gaussian mixture model (GMM). To our knowledge, this is the first result connecting a statistical model of counts to the distribution of cells in PCA space. Our derivation depends on methods from random matrix theory (RMT) and we extend and establish connections between RMT and scRNAseq count matrices. We then use the predicted GMM to test the metacell model on 11 datasets. The predictions are accurate for counts generated by the metacell model but underestimate the variance of real cells for most of our datasets. We show that this error largely arises from gene correlations within metacells, which the metacell model does not capture. Some of our datasets are fit relatively well by the metacell model, suggesting that the fit of the metacell model to true datasets depends on sufficiently dense sampling. Finally, we show that the strength of noise in the PCA space varies significantly across metacells, suggesting that statistical modeling of PCA embeddings can improve existing downstream workflows. Our results provide a framework for deriving and testing statistical models of PCA embeddings from statistical models of counts and identify a limitation of the metacell description for scRNAseq datasets.

Longitudinal non-Gaussian data are common in biomedical and microbiome studies, where responses may not be normally distributed. However, many existing dynamic graphical models are developed for Gaussian responses, while many non-Gaussian graphical models focus on cross-sectional dependence and do not explicitly estimate temporal transition effects. We propose a latent Gaussian dynamic graphical model for longitudinal non-Gaussian responses, including binary, count, and compositional outcomes. The proposed framework jointly estimates the contemporaneous conditional dependence network and the lag-one temporal dependence structure. We further develop a penalized Monte Carlo expectation maximization procedure for estimation. Simulation studies across multiple graph structures and outcome types show that the proposed method can recover both dependence structures, with more stable recovery for the lag-one temporal dependence structure and more structure-dependent recovery for the contemporaneous conditional dependence network. We apply the proposed method to longitudinal data from the Integrative Human Microbiome Project (HMP2) using a multinomial observation model. The results reveal persistent lag-one temporal dependence and contemporaneous conditional dependence among dominant microbial taxa, demonstrating the utility of the proposed framework for dynamic network analysis of longitudinal non-Gaussian biomedical data.

Graph embedding transforms complex networks into low-dimensional representations for analysis and downstream machine learning. Existing methods often rely on eigenvalue ordering, random-walk sampling, or neural optimization. We introduce MOSAIC, Multistep Operator Spectral Accumulation with Isolated Components, an unsupervised framework that constructs an accumulated diffusion operator from successive powers of a graph transition matrix. MOSAIC uses Krylov subspace analysis to identify dominant spectral components according to their convergence behavior. This approach emphasizes persistent and spectrally separated diffusion modes rather than selecting embedding coordinates solely by eigenvalue magnitude. We evaluated MOSAIC using a biomedical graph projection derived from the Data Distillery Knowledge Graph that contained 5,171 human genes and 10,147 phenotypes. Across embedding dimensions from 16 to 128, MOSAIC achieved node-type classification AUROCs from 0.9863 to 0.9926 and exceeded FastRP, Node2Vec, HashGNN, and GraphSAGE. For gene-phenotype link prediction, validation and test associations were withheld before embedding generation. Under progressive removal of all edges or specific edge classes, MOSAIC retained more node-label information than FastRP, including after extreme graph sparsification. In a focused DICER1 analysis, 106 of 137 documented phenotype associations received scores above 0.90, demonstrating concordance with known pleiotropic biology and prioritizing additional relationships for independent evaluation. These results support MOSAIC as a fully unsupervised framework for learning informative graph representations. Although demonstrated using a biomedical graph projection, MOSAIC is domain independent and can be applied to other complex networks

We present a lossless compression algorithm for indexing a protein database while supporting fast taxonomic classification in the method Centrifuger. The algorithm is a new scheme of the previously proposed run-block compression algorithm to reduce the size of the Ferragina-Manzini (FM) index, and it scales better with alphabet size than the original. On the RefSeq prokaryotic and viral protein sequences, Centrifuger reduces the memory footprint by over a third compared to the method Kaiju that builds on a plain FM-index, while having comparable running time. Furthermore, the compressed FM-index is lossless and can locate matches of arbitrary length, which helps Centrifuger achieve greater accuracy than Kraken2, a k-mer-based taxonomic classification method. We leverage the computational efficiency of Centrifuger to create an index of size 182 GB for classifying the reads against the full nr database that contains about 250 billion amino acid characters. Using this index, Centrifuger reveals different SARS-CoV-2 infection states and viral transcriptome profiles across human cell types from single-cell RNA-seq data.

Several tools are available to detect differential gene expression. The predominant approach of different statistical methods is to consider changes in gene expression independent from their genomic context. Previously, DELocal was introduced that takes into consideration gene expression dynamics in the genomic neighbourhood. Whereas DELocal was shown to identify developmentally regulated genes in developing mouse tooth, its performance in other biological conditions remains undetermined. To benchmark DElocal together with other methods, we develop another tool 'broadSeq'. It provides easy to use interface for methods like DELocal , deseq2, limma, edger, EBSeq and many more. Here we report benchmarking results on different data sets that comprise wild type developing mouse organs from different embryonic time points, and also human cancer data. The examined datasets are from the heart, kidney, liver, brain, forebrain, and hindbrain. We found that most methods produce relatively comparable lists of differentially expressed genes. Nevertheless, DELocal identifies partially distinct sets of genes compared to the other methods. Evaluating the identified genes using gene and disease ontology terms show that the uniquely detected genes by DELocal are biologically relevant and also highly expressed compared to the genes identified by the other methods. Because consistently highly expressed genes show relatively subtle fold changes, the neighbourhood concept of DELocal appears to correct for the difficulty in detecting these differential expression of highly expressed genes. Finally, we provide the tools as combined bioconductor packages, allowing researchers to explore and decide the best methods for their data and questions.

Microbial sequencing data are subject to sample-specific and species-specific distortions that complicate inference about true microbial abundances. Recent work established assumptions under which relative fold-differences in microbial abundances can be estimated from sequencing data, accompanied by an estimation and testing method called radEmu. While radEmu has excellent error rates, its computational cost can be prohibitive for large-scale analyses. Here we introduce fastEmu, a scalable approximation to radEmu's testing methodology. fastEmu uses a small, data-driven subset of taxa to approximate the computationally expensive test statistic. In a metagenomic analysis, we show a reduction in per-taxon runtime from 4-8 hours to 0.5-45 seconds, and simulations illustrate no adverse impact on estimation error or Type I error rate control. fastEmu enables testing differential abundance for thousands of microbial taxa while retaining clear parameter interpretation and error rate control.

Recent advances in artificial intelligence have enabled models to capture protein sequence and structural features. We present PeTriPPI2, a hybrid framework for protein--protein interaction (PPI) prediction that combines sequence embeddings from ESM-2 with structural representations from PeTriBERT, an encoder originally designed for inverse folding. We evaluate this sequence--structure approach on the Pinder dataset. On a test set of 2,342 protein pairs, the selected checkpoint achieves 0.898 accuracy, 0.917 precision, 0.876 recall, and 0.896 F1 score for the interacting class (AUROC 0.963). Compared with SpatialPPIv2 on the same test data, PeTriPPI2 has higher precision (0.917 versus 0.886) but lower recall (0.876 versus 0.942) and F1 score (0.896 versus 0.913). Component ablations indicate that both the sequence and structural pathways contribute to this configuration.

Biological populations have been evolving in multi-dimensional trait spaces. Ecological interactions can induce the evolutionary diversification of a biological population into distinct multiple populations, which is called evolutionary branching in adaptive dynamics theory. If possible mutations are restricted to certain directions due to developmental, physiological, or physical constraints, the resulting evolutionary dynamics may be restricted to subspaces with fewer dimensions than the original trait spaces. In general, such constraints are also evolving slowly, due to rare and slight mutational deviations from the constraint subspaces. In this sense, constraints are not strict. However, the method for describing non-strict constraints that is globally consistent has yet to be developed. This study proposes a general and flexible method for constructing a non-strict constraint subspace from the corresponding strict constraint subspace. This method enables efficient analysis of adaptive evolution including evolutionary branching under non-strict constraints in the framework of adaptive dynamics theory. Applying this method to a resource competition model shows that the shape of the carrying capacity distribution can be treated as a non-strict constraint curve, and that the likelihood of evolutionary branching is affected not only by the shape of the constraint curve but also by its slow evolution.

The advent of high-throughput sequencing technologies has made it essential to employ advanced tools for data integration and interpretation. In this work, we propose to build upon two existing efforts by expanding their scope and applicability to multi-omic data integration through the development of a pipeline that generates a topologically informed ranking of genes to assess gene relevance to disease. Specifically, we exploit topological data analysis, capitalizing on: its graph-based framework, to refine complex interactions characteristic of biological systems; its intrinsic explainability, to enable interpretation of model outputs in terms of biologically meaningful features; its robustness to small input perturbations, to ensure stable results in the presence of noisy data, as omics are. In particular, we employ Harmonic Persistent Homology (HPH), and propose two relevant advances. First, we adopt an integration and clustering tool (i.e. iNETgrate) as an informed means to reduce data size and integrate multiple omic layers by aggregating methylation beta values from loci to the gene level, to be combined with gene expression values. This enables the application of HPH, otherwise limited by the computational burden. Second, as a non-trivial consequence, we apply HPH to multi-omic molecular data, and not to patients as done so far, enlarging HPH scope. Finally, we follow up on recent efforts toward performance standardization, and validate our results against biomarkers from a well-known TCGA breast cancer benchmark dataset.

Biomarker discovery from high-dimensional RNA sequencing data remains challenging. Conventional methods such as differential expression analysis, pairwise protein-protein interaction networks, and weighted gene co-expression network analysis suffer from false-positive interactions, transitivity-driven noise, and arbitrary thresholds for feature selection. Addressing these limitations, this study proposed a novel framework, Hypergraph Regularised Fuzzy Rough Gene Selection (HR-FRGS), for knowledge-driven biomarker discovery. This pipeline integrates information-rich protein clustering with entropy filtering and PPI-pruned co-expression interactions to form a dual-layer hypergraph. A hypergraph-based scoring method that combines Random Walk with Restart diffusion and Hypergraph Betweenness Centrality is used to identify genes that are central in the biological network. Furthermore an autonomous fuzzy rough set selection eliminates arbitrary threshold dependency for gene selection. HR-FRGS was applied to four TCGA cancer cohorts: Lung Adenocarcinoma, Head and Neck Cancer, Kidney Clear Cell Carcinoma, and Colon Cancer, which reduced over 60,000 transcripts to compact biomarker panels. Validation using classical Machine Learning algorithms with cross-validation demonstrated that biomarker panels matched or exceeded the classification performance of the full transcriptome (AUC > 0.99, MCC > 0.94) and confirmed generalizability (MCC $\geq$ 0.893). Benchmarking against baseline methods and an ablation study showed the contribution of major parts of the algorithm. Functional enrichment analysis captures established pan-cancer hallmarks and cohort-specific oncogenic mechanisms, confirming the biological relevance of the identified genes.

In this work, we focus on predicting fold-switching phenomena observed in certain protein structures. Proteins that undergo fold-switching are called metamorphic proteins. We model fold-switch prediction in proteins as a position/residue-wise classification task. We use various features extracted from a protein's sequence as well as structure to train the classifiers. We propose a novel class of handcrafted parametrised features that encode the structural context of a residue. We also use embeddings from existing protein language models (PLMs) to train the classifiers. We train binary and 3-class classifiers, with the 3-class classifiers providing granular information on the type of fold-switch. In our computations, higher classification performance is achieved for binary prediction as compared to 3-class granular prediction. The handcrafted features perform comparably with the PLM embeddings on our standard test set. A rigorous error analysis by us during the review period revealed a critical data-leakage vulnerability in the metamorphic protein test data. This analysis also revealed that the PLM-based classifier suffers sequence-memorisation bias and fails to generalise to novel fold-switching sequences. On the other hand, one of our domain-informed handcrafted feature based classifier detects fold-switching without relying on sequence memory.

Abstract Accurate taxonomic classification in metagenomics is often compromised by low-complexity sequences, which lead to chance matches that in turn cause sequences to be misclassified. Here we present Kmask, an entropy-based masking tool implemented for use either standalone or as part of Kraken [1,2] database construction, which replaces low-entropy regions with Ns. Using a sliding window size aligned with Kraken's default k-mer length and parameters optimized across 12 control bacterial genomes spanning a broad range of GC content, Kmask efficiently removes low-complexity sequences while retaining high-complexity regions. To benchmark performance, we applied Kmask to a newly constructed database, Microbial2025, that contains over 71,000 bacterial, archaeal, viral, and fungal genomes, and we then classified human reads against both masked and unmasked versions of the database using KrakenUniq [2]. We found that Kmask substantially reduced misclassifications, driving down the false positive rate to 5.78% from 7.52%. Notably, Kmask performed comparably to an SDUST-masked [3] database, achieving a similar false positive rate (5.78% vs. 5.17%) while masking out fewer bases (1.33% vs. 1.85%). We also tested Kmask on a database of human cancer sequences, where we found that it eliminated many false positives caused by low-complexity matches between bacterial genomes and human DNA. These results demonstrate that Kmask is an effective method for masking low-complexity sequences in large microbial databases, thus improving the accuracy of metagenomic classification.

In network medicine, the protein-protein interaction network, or interactome, is an essential resource for identifying candidate proteins underlying diseases and other phenotypic traits. Indeed, it can be leveraged via the guilt-by-association paradigm, which asserts that interacting proteins are likely to participate in the same molecular processes. Network propagation that combines the topology of the interactome with prior knowledge about disease genes is a promising strategy to identify new candidate genes contributing to diseases. However, existing network propagation algorithms are often biased toward highly connected proteins, called "hubs". We present BFWalk, a novel network propagation algorithm designed to avoid inflated scores for hubs. We tested BFWalk across four human phenotypes: multiple morphological abnormalities of the sperm flagella, dyschromatopsia, hypertrophic cardiomyopathy, and chronic kidney disease. For each phenotype, we assessed whether BFWalk could recover known phenotype-associated genes and whether new candidate genes were enriched in relevant tissues. Depending on the phenotype, BFWalk outperformed or matched the state-of-the-art network propagation methods. Furthermore, the results confirmed that BFWalk is free of bias toward high-degree nodes. Therefore, BFWalk is a powerful method for identifying new genes across diverse phenotypes and constitutes a strong alternative to existing network propagation algorithms. Availability: BFWalk is implemented in Python and C. The code is available under the GNU GPL (v3.0) on: https://github.com/jedrzejkubica/BFWalk

Motivation: De novo genome assembly relies on accurate path selection through complex graphs, yet unitig graph branching points remain a major ambiguity source. Multiple extensions may be equally plausible from topology alone, but sequence-level consistency of competing paths remains largely unexplored. Model: We introduce RASTRO, a reference-free framework for scoring candidate paths at unitig branching points using nucleotide sequence information. RASTRO-S employs a sample-specific n-gram model to estimate extension plausibility via bits-per-base surprisal, while RASTRO-L uses the Evo2 pretrained DNA language model for autoregressive path scoring. Both operate directly on unitig graphs without reference genomes or assembled contigs. To evaluate branch decisions independently, we develop a read-based validation strategy using boundary k-mers spanning the junction between graph context and candidate extensions. Results: Experiments on HIV-1 and SARS-CoV-2 unitig graphs from Logan show sequence-based models provide meaningful match with read-supported choices. Among 20,000 sampled branching points, unique read-supported choices were identified for 15,319 HIV-1 branches (76.6%) and 19,882 SARS-CoV-2 branches (99.4%). RASTRO-S improved match from 54.1% to 75.8% for HIV-1 and from 58.6% to 76.4% for SARS-CoV-2 as context length rose from n=5 to n=10. RASTRO-L also improved with context length, reaching 49.96% and 49.25% match. Both outperformed random selection and the longest-extension greedy baseline. Conclusion: RASTRO provides an efficient, interpretable framework for evaluating ambiguous path decisions in unitig graphs using raw-read evidence rather than contig-level alignment. Results demonstrate the potential of lightweight, sample-adaptive sequence models for improving local decision-making in genome assembly workflows.